In [1]:
# Lee Mouse WT scRNA-seq data figures
# 05-20-2025
# Minhyung Kim
library(Seurat)
library(ggplot2)
library(dplyr)
set.seed(2025)
path <- "~/projects/2025/LeeNaeun/"
path_out <- "~/projects/2025/LeeNaeun/out/"
source("~/projects/MK_Rcode/volcano_scdata.R")
Loading required package: SeuratObject
Loading required package: sp
Attaching package: ‘SeuratObject’
The following objects are masked from ‘package:base’:
intersect, t
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
In [2]:
if (!file.exists(paste0(path,"Lee_WT_052025.rds"))) {
path_1 <- "~/projects/2025/LeeNaeun/HN00233422_result_10X/"
path_2 <- "~/projects/2025/LeeNaeun/HN00233483_result_10X/"
B1 <- c(1, 2)
B2 <- 3
f_name <- "/count/sample_filtered_feature_bc_matrix"
Lee_Mouse <- list()
CITE <- list()
for (i in B1) {
f_path <- paste0(path_1, i, "/per_sample_outs/", i, f_name)
data <- Read10X(data.dir = f_path)
Lee_Mouse[[paste0("sample_",i)]] <- CreateSeuratObject(counts = data$`Gene Expression`, project = paste0("sample_",i))
Lee_Mouse[[paste0("sample_",i)]]$batch <- "B1"
CITE[[paste0("sample_",i)]] <- data$`Antibody Capture`
}
for (i in B2) {
f_path <- paste0(path_2, i, "/per_sample_outs/", i, f_name)
data <- Read10X(data.dir = f_path)
Lee_Mouse[[paste0("sample_",i)]] <- CreateSeuratObject(counts = data$`Gene Expression`, project = paste0("sample_",i))
Lee_Mouse[[paste0("sample_",i)]]$batch <- "B2"
CITE[[paste0("sample_",i)]] <- data$`Antibody Capture`
}
Sobj <- merge(Lee_Mouse[["sample_1"]], y = list(Lee_Mouse[["sample_2"]], Lee_Mouse[["sample_3"]]), add.cell.ids = paste0("sample_",c(1:3)))
saveRDS(Sobj, file = paste0(path,"Lee_WT_052025.rds"))
saveRDS(CITE, file = paste0(path_out,"Lee_WT_CITE_052025.rds"))
rm(path_1, path_2, B1, B2, f_name, Lee_Mouse)
}
if (!exists("Sobj")) {
Sobj <- readRDS(paste0(path,"Lee_WT_052025.rds"))
CITE <- readRDS(paste0(path_out,"Lee_WT_CITE_052025.rds"))
}
In [3]:
CITE
[[ suppressing 34 column names ‘AAACCAAAGACAATTC-1’, ‘AAACCAAAGACGACGA-1’, ‘AAACCAAAGATGTATC-1’ ... ]] [[ suppressing 34 column names ‘AAACCAAAGACTATGT-1’, ‘AAACCAAAGCACAATT-1’, ‘AAACCAAAGCGAAGAT-1’ ... ]] [[ suppressing 34 column names ‘AAACCAAAGACTAATC-1’, ‘AAACCAAAGATTCGAT-1’, ‘AAACCAAAGCATCTTC-1’ ... ]]
$sample_1
7 x 20092 sparse Matrix of class "dgCMatrix"
CD138 1 . 2 3 1 1 3 4 1 1 1 1 4 2 2 10 4 2 7
CD19 109 177 124 114 187 69 129 162 21 22 13 160 83 235 135 148 324 151 77
CD11b 3 2 4 1 6 4 5 11 . 44 4 2 3 5 3 3 4 6 1
CD11c 8 4 5 1 4 6 3 5 4 16 1 8 4 2 2 3 9 23 3
CD43 14 17 23 14 12 84 38 50 22 105 4 8 25 56 54 59 257 18 28
CD5 26 16 20 18 26 49 31 38 18 18 16 12 18 30 84 20 263 24 12
CD21_35 39 133 88 42 58 28 52 89 8 1 1 191 42 63 137 68 12 55 21
CD138 1 1 6 2 5 13 2 9 3 2 4 1 1 6 5 ......
CD19 124 112 174 89 272 258 93 568 101 199 407 59 174 158 143 ......
CD11b 5 3 18 1 4 5 3 6 3 6 7 2 . 8 2 ......
CD11c 2 1 5 3 9 4 6 35 3 5 38 1 11 3 7 ......
CD43 33 79 23 26 63 59 7 65 16 29 58 10 13 53 23 ......
CD5 20 18 37 17 32 27 9 56 84 18 44 17 32 45 41 ......
CD21_35 63 77 43 31 273 190 34 808 48 179 1104 8 231 87 49 ......
.....suppressing 20058 columns in show(); maybe adjust options(max.print=, width=)
..............................
$sample_2
7 x 19107 sparse Matrix of class "dgCMatrix"
CD138 7 5 2 2 8 5 21 2 1 3 8 4 5 7 4 4 6 4
CD19 22 221 135 146 226 138 80 105 84 137 266 104 100 248 222 190 252 37
CD11b 58 558 3 5 . 4 . 3 1 1 3 4 2 8 3 3 3 1
CD11c 1194 12 1 5 2 2 3 3 1 4 1 . 1 17 12 3 4 2
CD43 75 595 36 8 6 33 33 33 5 39 154 24 14 111 43 206 20 .
CD5 244 73 28 17 23 37 23 54 15 23 49 24 25 202 47 87 38 18
CD21_35 11 77 54 124 181 72 49 75 46 35 94 8 77 57 124 4 306 67
CD138 1 7 12 3 7 3 5 9 2 1 1 4 1 11 5 4 ......
CD19 186 49 117 123 117 83 103 167 250 174 202 188 132 218 121 127 ......
CD11b 4 5 1 2 . 3 5 2 5 4 2 2 2 5 2 2 ......
CD11c 3 3 3 3 2 1 4 1 3 10 5 3 1 1 5 2 ......
CD43 9 18 21 41 10 4 34 21 23 27 15 35 13 21 19 23 ......
CD5 35 29 24 19 15 27 43 14 66 30 48 34 22 53 19 17 ......
CD21_35 159 17 34 59 130 44 48 62 86 67 238 179 101 176 90 105 ......
.....suppressing 19073 columns in show(); maybe adjust options(max.print=, width=)
..............................
$sample_3
7 x 24589 sparse Matrix of class "dgCMatrix"
CD138 1 6 21 13 7 6 2 8 13 6 7 97 5 6 4 . 2 15
CD19 84 152 221 569 176 81 124 183 113 161 100 150 134 14 201 131 155 189
CD11b 7 7 5 133 6 2 6 44 13 5 3 2 8 1 6 2 6 10
CD11c 1 2 5 9 3 1 3 5 2 1 5 3 4 2 7 3 2 7
CD43 13 36 55 522 29 5 15 70 10 14 53 31 19 7 13 4 6 32
CD5 18 14 30 245 30 109 19 40 27 18 94 18 15 13 24 11 43 44
CD21_35 2 125 31 39 122 17 75 106 59 73 15 73 69 4 60 73 367 99
CD138 6 4 7 12 7 3 6 11 2 7 7 12 3 3 13 7 ......
CD19 161 99 137 27 33 16 111 344 105 145 15 254 196 111 11 39 ......
CD11b 13 3 4 294 1926 7 4 15 2 4 32 9 3 2 5 9 ......
CD11c 6 2 3 257 67 2 8 3 3 4 7 2 3 . 5 5 ......
CD43 7 25 8 617 3304 5 30 165 6 10 120 75 1 19 9 13 ......
CD5 81 11 16 52 56 11 26 53 26 15 16 26 17 19 6 28 ......
CD21_35 322 86 38 14 12 7 24 112 32 67 4 201 133 45 12 21 ......
.....suppressing 24555 columns in show(); maybe adjust options(max.print=, width=)
..............................
In [4]:
Sobj <- JoinLayers(Sobj)
dim(CITE$sample_1)
dim(CITE$sample_2)
dim(CITE$sample_3)
- 7
- 20092
- 7
- 19107
- 7
- 24589
In [5]:
bb <- cbind(
CITE$sample_1,
CITE$sample_2,
CITE$sample_3
)
colnames(bb) <- colnames(Sobj)
In [6]:
length(colnames(Sobj))
head(colnames(Sobj))
63788
- 'sample_1_AAACCAAAGACAATTC-1'
- 'sample_1_AAACCAAAGACGACGA-1'
- 'sample_1_AAACCAAAGATGTATC-1'
- 'sample_1_AAACCAAAGCCAACCG-1'
- 'sample_1_AAACCAAAGCCGATTA-1'
- 'sample_1_AAACCAAAGCTGCGTG-1'
In [7]:
Sobj
An object of class Seurat 33696 features across 63788 samples within 1 assay Active assay: RNA (33696 features, 0 variable features) 1 layer present: counts
In [8]:
# QC
Sobj[["percent.mt"]] <- PercentageFeatureSet(Sobj, pattern = "^mt-")
In [9]:
VlnPlot(Sobj, features = c("nFeature_RNA"), pt.size = 0.1)
VlnPlot(Sobj, features = c("nFeature_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("nCount_RNA"), pt.size = 0.1)
VlnPlot(Sobj, features = c("nCount_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("percent.mt"), pt.size = 0.1)
VlnPlot(Sobj, features = c("percent.mt"), pt.size = 0)
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.” Warning message: “The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0. ℹ Please use the `layer` argument instead. ℹ The deprecated feature was likely used in the Seurat package. Please report the issue at <https://github.com/satijalab/seurat/issues>.” Warning message: “`PackageCheck()` was deprecated in SeuratObject 5.0.0. ℹ Please use `rlang::check_installed()` instead. ℹ The deprecated feature was likely used in the Seurat package. Please report the issue at <https://github.com/satijalab/seurat/issues>.” Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
In [10]:
p1 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "percent.mt", raster = F)
p2 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F)
p1+p2
t_1 <- table(Sobj$orig.ident)
t_1
sample_1 sample_2 sample_3 20092 19107 24589
In [11]:
cutoffs <- data.frame(matrix(NA,nrow = 3, ncol = 3))
colnames(cutoffs) <- c("nF_cutoff", "mt_cutoff", "mt_cutoff_sd")
rownames(cutoffs) <- paste0("sample_", c(1:3))
cell_ind <- list()
for (a in rownames(cutoffs)) {
cutoffs[a,"nF_cutoff"] <- median(Sobj$nFeature_RNA[Sobj$orig.ident==a]) + sd(Sobj$nFeature_RNA[Sobj$orig.ident==a])*5
cutoffs[a,"mt_cutoff_sd"] <- median(Sobj$percent.mt[Sobj$orig.ident==a]) + sd(Sobj$percent.mt[Sobj$orig.ident==a])*2
}
In [12]:
cutoffs[,"mt_cutoff"] <- cutoffs[,"mt_cutoff_sd"]
cutoffs[cutoffs[,"mt_cutoff_sd"]>25,"mt_cutoff"] <- 25
cell_ind <- NULL
for (s_name in rownames(cutoffs)) {
tm_ind <- colnames(Sobj[,Sobj$orig.ident==s_name & Sobj$nFeature_RNA < cutoffs[s_name,"nF_cutoff"] & Sobj$nFeature_RNA>300 & Sobj$percent.mt < cutoffs[s_name,"mt_cutoff"]])
cell_ind <- c(cell_ind,tm_ind)
}
In [13]:
expressed_cell_count <- data.frame(matrix(NA,nrow = dim(Sobj)[1], ncol = 9))
colnames(expressed_cell_count) <- rownames(cutoffs)
for (i in rownames(cutoffs)) {
expressed_cell_count[,i] <- rowSums(Sobj@assays$RNA$counts[,Sobj$orig.ident==i]>1)
}
print(table(Sobj$orig.ident))
sample_1 sample_2 sample_3 20092 19107 24589
In [14]:
print(colSums(expressed_cell_count>10))
keep_genes <- rowSums(expressed_cell_count>10)>0
sample_1 sample_2 sample_3 <NA> <NA> <NA> <NA> <NA>
12635 12425 12374 NA NA NA NA NA
<NA>
NA
In [15]:
Sobj$cell_filtered <- is.element(colnames(Sobj), cell_ind)
for (s_name in rownames(cutoffs)) {
a0 <- VlnPlot(Sobj[,Sobj$orig.ident==s_name], features = c("nFeature_RNA","nCount_RNA","percent.mt"), pt.size = 0)
a1 <- FeatureScatter(Sobj[,Sobj$orig.ident==s_name], feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "cell_filtered") + theme(legend.position = "none") + ggtitle(s_name)
a2<- FeatureScatter(Sobj[,Sobj$orig.ident==s_name], feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "cell_filtered") + theme(legend.position = "none") + ggtitle(s_name)
a3<- FeatureScatter(Sobj[,Sobj$orig.ident==s_name], feature1 = "nFeature_RNA", feature2 = "percent.mt", group.by = "cell_filtered") + theme(legend.position = "none") + ggtitle(s_name)
print(a0)
print(a1+a2+a3)
}
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
In [16]:
Sobj <- subset(Sobj, subset = cell_filtered, features = rownames(Sobj)[keep_genes])
dim(Sobj)
- 33696
- 62860
In [17]:
Sobj[["CITE"]] <- CreateAssayObject(counts = bb[,cell_ind])
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
Warning message:
“Feature names cannot have underscores ('_'), replacing with dashes ('-')”
In [18]:
t_2 <- table(Sobj$orig.ident)
t_2
sample_1 sample_2 sample_3 19869 18861 24130
In [19]:
VlnPlot(Sobj, features = c("nFeature_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("nCount_RNA"), pt.size = 0)
VlnPlot(Sobj, features = c("percent.mt"), pt.size = 0)
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.” Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
Warning message: “Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.”
In [20]:
Idents(Sobj) <- Sobj$orig.ident
p1 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "percent.mt", raster = F)
p2 <- FeatureScatter(Sobj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F)
p1+p2
In [21]:
Sobj <- NormalizeData(Sobj, normalization.method = "LogNormalize", scale.factor = 10000)
Sobj <- FindVariableFeatures(Sobj, selection.method = "vst", nfeatures = 2000)
Sobj <- ScaleData(Sobj)
Sobj <- RunPCA(Sobj)
DimPlot(Sobj)
ElbowPlot(Sobj)
Normalizing layer: counts Finding variable features for layer counts Centering and scaling data matrix PC_ 1 Positive: Igkc, Iglc3, Ighm, Gm30211, Cr2, Mzb1, Kcnq5, Myc, Mast4, E330020D12Rik Immp2l, Inpp4b, Maml2, 2010309G21Rik, Iglc1, Ncl, Hs3st1, Dtx1, Arhgap24, Rab30 Plk2, Cacna1e, Txndc5, Mif, Ly6a, Hdac9, Ift140, Rps2, Pik3r4, Iglv3 Negative: Fcer1g, Ifitm3, Lyz2, Cd300a, Ifitm2, Csf1r, Alox5ap, Ccl6, Fyb, Clec7a Ifitm6, S100a4, Metrnl, Clec4a1, Clec4a3, Igsf6, Lgals3, Emilin2, Tnfrsf1a, Gm36161 Sirpa, Gpr141, Cd300c2, Pid1, Gsr, Lrp1, Gm21188, Ifi207, Itgam, Lst1 PC_ 2 Positive: Ace, Pparg, Treml4, Ear2, Cd300e, Cd300ld, Adgre4, S1pr5, Fabp4, Eno3 Pilra, Tgm2, Smpdl3b, Fcgr4, Zfyve9, Slc11a1, Apoc2, Rap1gap2, Pilrb2, Tiam2 Hfe, Arhgef10l, Cd300c2, Trem3, Gcnt2, Havcr2, Pglyrp1, Fgd4, Adrb1, Lilra5 Negative: Mki67, Pclaf, Birc5, Kif11, Ccna2, Top2a, Cdca8, Cdk1, Spc24, Tpx2 Cks1b, Cdca3, Racgap1, Stmn1, Kif15, Ube2c, Hmmr, Cenpf, Ccnb2, Prc1 Kif22, H2ac4, Ncapg, H2az1, Clspn, Esco2, Neil3, Bub1, Kif4, Cdca5 PC_ 3 Positive: F13a1, Fn1, Ccr2, AI839979, Ccl9, Cd177, Ly6c2, Stxbp6, Plcb1, Slfn1 C3, Hp, Ly6a2, Wfdc17, Clec4a2, Gas7, Olfm1, Hacd4, Ifi211, Ms4a8a Msr1, Cd14, S100a6, Thbs1, Mmp8, Mcub, Tgfbi, Mgst1, Arhgef37, S100a4 Negative: Pparg, Cd300e, Ace, Treml4, Eno3, Ear2, Tgm2, Cd9, Lair1, Cd300ld Havcr2, Cd36, S1pr5, Fabp4, Adgre4, Gcnt2, Itgax, Fcgr4, Spn, Hfe Cdk14, Mki67, Grk3, Pilra, Smpdl3b, Slc12a2, Pclaf, Rap1gap2, Runx2, Ppm1h PC_ 4 Positive: Mki67, Pclaf, Birc5, Ccna2, Kif11, Cdk1, Spc24, Cdca3, Top2a, Cdca8 Tpx2, Hmmr, Ube2c, Cenpf, Kif15, Ccnb2, Prc1, Kif22, H2ac4, Stmn1 Esco2, Cenpe, Neil3, Clspn, Bub1, Nuf2, Kif4, Rrm2, Ncapg, Cks1b Negative: Cd7, Klrd1, Tcf7, Itk, Skap1, Nkg7, Ms4a4b, Cd3d, Cd3g, Gm2682 Cd3e, Lat, Cd28, Prkcq, Trbc2, Cd247, Il2rb, Ccl5, Cd27, Slc9a9 Il7r, Ccnd1, Klrk1, Srgap3, Prkch, Thy1, Ctsw, Tbc1d4, Tox, Bcl11b PC_ 5 Positive: Atxn1, Pik3r4, Zc3h12c, Ahr, Myof, Fcrl5, Dtx1, Cr2, Myc, Atf3 Nebl, S1pr3, Plac8, Tank, Ackr3, Odc1, Dph5, Hs3st1, Cd86, Pdia4 Rilpl2, Tbc1d9, Camkmt, AW011738, Dnase1l3, Immp2l, Cybb, Hpgds, Marcks, Jun Negative: Itk, Tcf7, Ms4a4b, Nkg7, Cd3d, Gm2682, Cd3e, Cd3g, Trbc2, Lat Cd247, Prkcq, Cd28, Cd27, Skap1, Tox, Prkch, Il2rb, Lef1, Thy1 Bcl11b, Trbc1, Themis, Ctsw, Il7r, Camk4, Sh2d1a, Ctla2a, Txk, Ikzf2
In [22]:
Sobj <- RunUMAP(Sobj, dims = 1:15)
Warning message: “The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation' This message will be shown once per session” 17:01:14 UMAP embedding parameters a = 0.9922 b = 1.112 17:01:14 Read 62860 rows and found 15 numeric columns 17:01:14 Using Annoy for neighbor search, n_neighbors = 30 17:01:14 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 17:01:22 Writing NN index file to temp file /tmp/RtmpiDhTyy/file5b1f448442d2a 17:01:22 Searching Annoy index using 1 thread, search_k = 3000 17:01:46 Annoy recall = 100% 17:01:47 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 17:01:49 Initializing from normalized Laplacian + noise (using RSpectra) 17:01:55 Commencing optimization for 200 epochs, with 2692262 positive edges 17:01:55 Using rng type: pcg 17:02:36 Optimization finished
In [23]:
DimPlot(Sobj, raster = F)
In [24]:
library(harmony)
Loading required package: Rcpp
In [25]:
Sobj <- RunHarmony(Sobj,"orig.ident", plot_convergence = T)
harmony_embeddings <- Embeddings(Sobj, "harmony")
harmony_embeddings[1:5, 1:5]
Transposing data matrix Initializing state using k-means centroids initialization Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Warning message: “Quick-TRANSfer stage steps exceeded maximum (= 3143000)” Harmony 1/10 Harmony 2/10 Harmony 3/10 Harmony converged after 3 iterations Warning message: “The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0. ℹ Please use the `layer` argument instead. ℹ The deprecated feature was likely used in the Seurat package. Please report the issue at <https://github.com/satijalab/seurat/issues>.”
| harmony_1 | harmony_2 | harmony_3 | harmony_4 | harmony_5 | |
|---|---|---|---|---|---|
| sample_1_AAACCAAAGACAATTC-1 | 1.449980 | -1.9607286 | -2.1336740 | -0.1322129 | 3.99104252 |
| sample_1_AAACCAAAGACGACGA-1 | 2.540081 | 0.9835546 | 1.1427124 | 0.3662081 | -0.06780086 |
| sample_1_AAACCAAAGATGTATC-1 | 2.281648 | 1.0801510 | -0.4316720 | -1.0081381 | 2.18528764 |
| sample_1_AAACCAAAGCCAACCG-1 | 3.344910 | 1.1810851 | 0.1232371 | 1.0964129 | -0.49868922 |
| sample_1_AAACCAAAGCCGATTA-1 | 1.612576 | 0.1193344 | 0.3849787 | -0.7293387 | 1.97029635 |
In [26]:
DimPlot(Sobj, reduction = "harmony", group.by = "orig.ident", raster = F)
In [27]:
Sobj <- RunUMAP(Sobj, reduction = "harmony", dims = 1:15)
Sobj <- FindNeighbors(Sobj, reduction = "harmony", dims = 1:15)
Sobj <- FindClusters(Sobj, resolution = 0.5)
17:03:56 UMAP embedding parameters a = 0.9922 b = 1.112 17:03:56 Read 62860 rows and found 15 numeric columns 17:03:56 Using Annoy for neighbor search, n_neighbors = 30 17:03:56 Building Annoy index with metric = cosine, n_trees = 50 0% 10 20 30 40 50 60 70 80 90 100% [----|----|----|----|----|----|----|----|----|----| * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * | 17:04:03 Writing NN index file to temp file /tmp/RtmpiDhTyy/file5b1f4286fb1ac 17:04:03 Searching Annoy index using 1 thread, search_k = 3000 17:04:28 Annoy recall = 100% 17:04:29 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30 17:04:31 Initializing from normalized Laplacian + noise (using RSpectra) 17:04:36 Commencing optimization for 200 epochs, with 2705776 positive edges 17:04:36 Using rng type: pcg 17:05:17 Optimization finished Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 62860 Number of edges: 1625800 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.8656 Number of communities: 17 Elapsed time: 21 seconds
In [28]:
DimPlot(Sobj, raster = F)
DimPlot(Sobj, raster = F, group.by = "orig.ident")
In [29]:
Sobj <- ScaleData(Sobj, features = rownames(Sobj))
Centering and scaling data matrix Warning message: “Different features in new layer data than already exists for scale.data”
In [30]:
source("~/projects/MK_Rcode/setZ.R")
ABC_marker_mouse <- c("Cd19", "Cd80", "Cd86", "Fas", "Fcrla", "Fcrlb", "Fgl2", "Cxcl10", "Cxcr3", "Nkg7", "Itgam", "Itgb2l", "Itgb2", "Tbx21", "Zbtb32", "Fcer2a", "Cr2")
length(ABC_marker_mouse)
is.element(ABC_marker_mouse,rownames(Sobj))
ABC_marker_mouse <- ABC_marker_mouse[is.element(ABC_marker_mouse,rownames(Sobj))]
length(ABC_marker_mouse)
17
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
- TRUE
17
In [31]:
Sobj[["ABC_score"]] <- setZ(ABC_marker_mouse,rownames(Sobj@assays$RNA$scale.data),Sobj@assays$RNA$scale.data)
In [32]:
LSP1_gene_mouse <- c("Cebpa", "Cebpb", "Il1b", "Ccl3", "Tnf", "Ccl4", "Ccl5", "Spp1", "Cd7", "Msr1", "Itgam", "Kit", "Cd14", "Csf3r", "Il1r2", "Csf1r", "Itga1", "Itga5", "Itgb5", "Itgb1", "Camp", "Thbs1", "C3", "Mefv", "Fos", "Cybb", "Card9", "Irf7", "Oas3", "Nlrp3", "Ctsb", "Ikbke", "Mpo")
length(LSP1_gene_mouse)
33
In [33]:
LSP1_gene_sub <- LSP1_gene_mouse[is.element(LSP1_gene_mouse,rownames(Sobj))]
length(LSP1_gene_sub)
33
In [34]:
Sobj[["LSP1_score"]] <- setZ(LSP1_gene_sub,rownames(Sobj@assays$RNA$scale.data),Sobj@assays$RNA$scale.data)
In [35]:
Sobj$LSP1_score_pos <- Sobj$LSP1_score > 0
Sobj$ABC_score_pos <- Sobj$ABC_score > 0
FeaturePlot(Sobj, features = "LSP1_score") + ggtitle("LSP1 score")
FeaturePlot(Sobj, features = "ABC_score") + ggtitle("ABC score")
p0 <- DimPlot(Sobj, group.by = "LSP1_score_pos") + ggtitle("LRB")
p0
p0_1 <- DimPlot(Sobj, group.by = "ABC_score_pos") + ggtitle("ABC")
p0_1
In [36]:
if (!file.exists(paste0(path_out,"Lee_WT_harmony_052025.rds"))) {
saveRDS(Sobj, paste0(path_out,"Lee_WT_harmony_052025.rds"))
saveRDS(Sobj@meta.data, paste0(path_out,"Lee_WT_metadata_052025.rds"))
}
In [37]:
marker_list <- c("Myc", "Ebf1", "Pax5", "Dtx1", "Bank1", "Cd55", "Cd19", "Cd93", "Cd24a", "Cd38", "Ighm", "Ighd", "Ighg", "Cd23", "Cd21", "Sell", "Cd5", "Cd43", "Cd27", "Cd80", "Cd73", "Pdcd1lg2", "Bcl6", "Aicda", "Mki67", "Cxcr4", "Cd83", "Cd86", "Cd74", "Prdm1", "Xbp1", "Irf4", "Cd138", "Jchain", "Mzb1")
In [38]:
Idents(Sobj) <- Sobj$seurat_clusters
pdf(paste0(path_out,"marker_dotplot_052125.pdf"), width = 20, height = 6)
DotPlot(Sobj, features = marker_list) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)
dev.off()
Warning message: “The following requested variables were not found: Ighg, Cd23, Cd21, Cd43, Cd73, Cd138” Scale for colour is already present. Adding another scale for colour, which will replace the existing scale.
pdf: 2
In [39]:
if (!file.exists(paste0(path_out,"marker_list_WT_052125.rds"))) {
all.markers <- FindAllMarkers(Sobj, test.use = "MAST", min.pct = 0.1, logfc.threshold = 0.58, verbose = F, only.pos = T)
saveRDS(all.markers, file = paste0(path_out,"marker_list_WT_052125.rds"))
}
if (!exists("all.markers")) {
all.markers <- readRDS(paste0(path_out,"marker_list_WT_052125.rds"))
}
In [40]:
marker_200 <- all.markers %>%
group_by(cluster) %>%
filter(avg_log2FC > 0.58) %>%
slice_head(n = 200) %>%
ungroup()
In [41]:
write.csv(marker_200, paste0(path_out, "top200_marker_genes.csv"))
In [42]:
Idents(Sobj) <- Sobj$seurat_clusters
Sobj <- RenameIdents(Sobj,
"0" = "Naive Follicular B", "1" = "Naive Follicular B", "2" = "Naive Follicular B", "3" = "Naive Follicular B",
"4" = "Transitional B (T2)", "5" = "Naive Follicular B", "6" = "Late transitional B", "7" = "Naive Follicular B",
"8" = "Plasma cell", "9" = "Myeloid-like GC B", "10" = "Myeloid-like GC B", "11" = "Myeloid-like GC B",
"12" = "Mki67+ GC DZ B", "13" = "T/B Doublet", "14" = "GC LZ B", "15" = "Naive Follicular B", "16" = "Plasmablast"
)
Sobj$annotation <- Idents(Sobj)
In [43]:
d1 <- DimPlot(Sobj, group.by = "seurat_clusters", label = T, repel = T)
d2 <- DimPlot(Sobj, group.by = "annotation", label = T, repel = T)
d1
d2
pdf(paste0(path_out,"umap_052125.pdf"), width = 6, height = 4)
d1 + ggtitle("")
d2 + ggtitle("")
dev.off()
pdf(paste0(path_out,"umap_053025v1.pdf"), width = 4, height = 4)
d1 + ggtitle("") + NoLegend()
d2 + ggtitle("") + NoLegend()
DimPlot(Sobj, group.by = "annotation", label = F) + ggtitle("") + NoLegend()
p0
p0 + ggtitle("") + NoLegend()
p0_1
p0_1 + ggtitle("") + NoLegend()
dev.off()
pdf: 2
pdf: 2
In [44]:
RidgePlot(Sobj, features = c("LSP1_score", "ABC_score"), ncol = 2)
d3 <- RidgePlot(Sobj, features = c("LSP1_score")) + NoLegend() + ggtitle("")
d4 <- RidgePlot(Sobj, features = c("ABC_score")) + NoLegend() + ggtitle("")
Picking joint bandwidth of 0.385 Picking joint bandwidth of 0.355
In [45]:
pdf(paste0(path_out, "RidgePlot_LSP1_ABC_052225.pdf"), width = 4, height = 4)
d3
d4
dev.off()
Picking joint bandwidth of 0.385 Picking joint bandwidth of 0.355
pdf: 2
In [46]:
celltype.prop <- as.data.frame(table(Sobj$annotation))
celltype.prop <- celltype.prop %>%
rename(celltype = Var1) %>%
mutate(Percent = Freq / sum(Freq)*100)
d5 <- ggplot(celltype.prop, aes(x = "", y = Percent, fill = celltype))+
geom_bar(stat = "identity", width = 0.9) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), axis.title.x = element_blank(), plot.margin = margin(10, 10, 20, 10)) +
NoLegend()
d5
In [47]:
pdf(paste0(path_out,"stacked_bar_WT_052225.pdf"), width = 1.5, height = 4)
d5
dev.off()
pdf: 2
In [48]:
celltype.prop
| celltype | Freq | Percent |
|---|---|---|
| <fct> | <int> | <dbl> |
| Naive Follicular B | 49063 | 78.0512249 |
| Transitional B (T2) | 5054 | 8.0400891 |
| Late transitional B | 2659 | 4.2300350 |
| Plasma cell | 1801 | 2.8650970 |
| Myeloid-like GC B | 2930 | 4.6611518 |
| Mki67+ GC DZ B | 488 | 0.7763283 |
| T/B Doublet | 413 | 0.6570156 |
| GC LZ B | 342 | 0.5440662 |
| Plasmablast | 110 | 0.1749920 |
In [49]:
table(Sobj$LSP1_score_pos)
table(Sobj$ABC_score_pos)
FALSE TRUE 47330 15530
FALSE TRUE 34780 28080
In [50]:
LSP1_pos_ind <- Sobj$LSP1_score_pos==TRUE & Sobj$ABC_score_pos==FALSE
ABC_pos_ind <- Sobj$LSP1_score_pos==FALSE & Sobj$ABC_score_pos==TRUE
DP_ind <- Sobj$LSP1_score_pos==TRUE & Sobj$ABC_score_pos==TRUE
DN_ind <- Sobj$LSP1_score_pos==FALSE & Sobj$ABC_score_pos==FALSE
Sobj$anno_ABC_LSP1 <- NA
Sobj$anno_ABC_LSP1[LSP1_pos_ind] <- "LRB"
Sobj$anno_ABC_LSP1[ABC_pos_ind] <- "ABC"
Sobj$anno_ABC_LSP1[DP_ind] <- "LRB & ABC"
Sobj$anno_ABC_LSP1[DN_ind] <- "Classical B"
unique(Sobj$anno_ABC_LSP1)
Idents(Sobj) <- factor(Sobj$anno_ABC_LSP1, levels = c("LRB", "ABC", "LRB & ABC", "Classical B"))
Sobj$anno_ABC_LSP1 <- Idents(Sobj)
levels(Idents(Sobj))
- 'Classical B'
- 'ABC'
- 'LRB'
- 'LRB & ABC'
- 'LRB'
- 'ABC'
- 'LRB & ABC'
- 'Classical B'
In [51]:
options(repr.plot.width=10, repr.plot.height=8)
Idents(Sobj) <- Sobj$anno_ABC_LSP1
levels(Sobj$anno_ABC_LSP1)
Idents(Sobj) <- Sobj$anno_ABC_LSP1
my_colors <- c("LRB" = "#ff0000", "ABC" = "#ffff00", "LRB & ABC" = "#00b0ff", "Classical B" = "#BFBFBF")
d6 <- DimPlot(Sobj, reduction = "umap", cols = my_colors)
d6
- 'LRB'
- 'ABC'
- 'LRB & ABC'
- 'Classical B'
In [52]:
if (!file.exists(paste0(path_out,"Lee_WT_harmony_080425.rds"))) {
saveRDS(Sobj, paste0(path_out,"Lee_WT_harmony_080425.rds"))
}
In [53]:
library(ggpointdensity)
library(scales)
In [54]:
df <- FetchData(Sobj, vars = c("umap_1", "umap_2", "LSP1_score_pos", "ABC_score_pos", "anno_ABC_LSP1"))
df_LSP1 <- df %>% filter(LSP1_score_pos)
df_ABC <- df %>% filter(ABC_score_pos)
In [55]:
options(repr.plot.width=6, repr.plot.height=5)
t1 <- ggplot(df_LSP1, aes(x = umap_1, y = umap_2)) +
geom_pointdensity(method = "neighbors", size = 0.1) +
scale_color_gradientn(
colors = c("grey", "grey", "red"),
values = rescale(c(0, 350, 1200)),
limits = c(0, 1200)
) +
theme_classic()
t2 <- ggplot(df_ABC, aes(x = umap_1, y = umap_2)) +
geom_pointdensity(method = "neighbors", size = 0.1) +
scale_color_gradientn(
colors = c("grey", "grey", "red"),
values = rescale(c(0, 350, 1200)),
limits = c(0, 1200)
) +
theme_classic()
t1
t2
In [56]:
pdf(paste0(path_out,"LRB_ABC_umap_density_080425.pdf"), width = 6, height = 5)
t1
t2
dev.off()
pdf: 2
In [57]:
t1_1 <- t1 +
labs(x = "", y = "") +
theme(
legend.position = "none",
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_line(),
axis.ticks.y = element_line()
)
t2_1 <- t2 +
labs(x = "", y = "") +
theme(
legend.position = "none",
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_line(),
axis.ticks.y = element_line()
)
t1_1
t2_1
In [70]:
tiff(filename = paste0(path_out,"LRB_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t1_1)
dev.off()
tiff(filename = paste0(path_out,"ABC_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t2_1)
dev.off()
pdf: 2
pdf: 2
In [69]:
path_out
'~/projects/2025/LeeNaeun/out/'
In [59]:
library(clusterProfiler)
library(enrichplot)
library(msigdbr)
library(stringr)
set.seed(2025)
hallmark_geneset <- msigdbr(species = "mouse", category = "H")
clusterProfiler v4.14.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
Please cite:
G Yu. Thirteen years of clusterProfiler. The Innovation. 2024,
5(6):100722
Attaching package: ‘clusterProfiler’
The following object is masked from ‘package:stats’:
filter
enrichplot v1.26.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
Please cite:
Guangchuang Yu. Using meshes for MeSH term enrichment and semantic
analyses. Bioinformatics. 2018, 34(21):3766-3767,
doi:10.1093/bioinformatics/bty410
Using human MSigDB with ortholog mapping to mouse. Use `db_species = "MM"` for mouse-native gene sets.
This message is displayed once per session.
Warning message:
“The `category` argument of `msigdbr()` is deprecated as of msigdbr 10.0.0.
ℹ Please use the `collection` argument instead.”
In [60]:
hallmark_geneset$gs_name_short <- str_extract(hallmark_geneset$gs_name, "(?<=HALLMARK_).*$")
head(hallmark_geneset)
| gene_symbol | ncbi_gene | ensembl_gene | db_gene_symbol | db_ncbi_gene | db_ensembl_gene | source_gene | gs_id | gs_name | gs_collection | ⋯ | gs_url | db_version | db_target_species | ortholog_taxon_id | ortholog_sources | num_ortholog_sources | entrez_gene | gs_cat | gs_subcat | gs_name_short |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | ⋯ | <chr> | <chr> | <chr> | <int> | <chr> | <dbl> | <chr> | <chr> | <chr> | <chr> |
| Abca1 | 11303 | ENSMUSG00000015243 | ABCA1 | 19 | ENSG00000165029 | ABCA1 | M5905 | HALLMARK_ADIPOGENESIS | H | ⋯ | 2025.1.Hs | HS | 10090 | EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|Treefam | 11 | 11303 | H | ADIPOGENESIS | ||
| Abcb8 | 74610 | ENSMUSG00000028973 | ABCB8 | 11194 | ENSG00000197150 | ABCB8 | M5905 | HALLMARK_ADIPOGENESIS | H | ⋯ | 2025.1.Hs | HS | 10090 | EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|PhylomeDB|Treefam | 12 | 74610 | H | ADIPOGENESIS | ||
| Acaa2 | 52538 | ENSMUSG00000036880 | ACAA2 | 10449 | ENSG00000167315 | ACAA2 | M5905 | HALLMARK_ADIPOGENESIS | H | ⋯ | 2025.1.Hs | HS | 10090 | EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|Treefam | 11 | 52538 | H | ADIPOGENESIS | ||
| Acadl | 11363 | ENSMUSG00000026003 | ACADL | 33 | ENSG00000115361 | ACADL | M5905 | HALLMARK_ADIPOGENESIS | H | ⋯ | 2025.1.Hs | HS | 10090 | EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|PhylomeDB|Treefam | 12 | 11363 | H | ADIPOGENESIS | ||
| Acadm | 11364 | ENSMUSG00000062908 | ACADM | 34 | ENSG00000117054 | ACADM | M5905 | HALLMARK_ADIPOGENESIS | H | ⋯ | 2025.1.Hs | HS | 10090 | EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|PhylomeDB|Treefam | 12 | 11364 | H | ADIPOGENESIS | ||
| Acads | 11409 | ENSMUSG00000029545 | ACADS | 35 | ENSG00000122971 | ACADS | M5905 | HALLMARK_ADIPOGENESIS | H | ⋯ | 2025.1.Hs | HS | 10090 | EggNOG|Ensembl|HGNC|HomoloGene|Inparanoid|NCBI|OMA|OrthoDB|OrthoMCL|Panther|Treefam | 11 | 11409 | H | ADIPOGENESIS |
In [63]:
if (!file.exists(paste0(path_out,"LSP1_vs_ABC.marker.rds"))) {
LSP1_vs_ABC.marker <- FindMarkers(Sobj, ident.1 = "LRB", ident.2 = "ABC", min.pct = 0.1, logfc.threshold = 0.1, test.use = "MAST")
saveRDS(LSP1_vs_ABC.marker, file = paste0(path_out,"LSP1_vs_ABC.marker.rds"))
}
if (!exists("LSP1_vs_ABC.marker")) {
LSP1_vs_ABC.marker <- readRDS(paste0(path_out,"LSP1_vs_ABC.marker.rds"))
}
In [64]:
sort_glist <- LSP1_vs_ABC.marker$avg_log2FC
names(sort_glist) <- rownames(LSP1_vs_ABC.marker)
sort_glist <- sort(sort_glist, decreasing = T)
gsea_result.LSP1vsABS <- GSEA(geneList = sort_glist, TERM2GENE = select(hallmark_geneset, gs_name_short, gene_symbol), seed = T, pvalueCutoff = 0.99)
using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019). preparing geneSet collections... GSEA analysis... leading edge analysis... done...
In [65]:
aa <- dotplot(gsea_result.LSP1vsABS, showCategory=10, split=".sign") + facet_grid(.~.sign)
aa
saveRDS(gsea_result.LSP1vsABS, file = paste0(path,"gsea_.rds"))
In [66]:
head(gsea_result.LSP1vsABS@result)
write.csv(gsea_result.LSP1vsABS@result, file = paste0(path_out,"gsea_mouse_LRB_vs_ABS.csv"))
| ID | Description | setSize | enrichmentScore | NES | pvalue | p.adjust | qvalue | rank | leading_edge | core_enrichment | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <int> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | <chr> | |
| INTERFERON_GAMMA_RESPONSE | INTERFERON_GAMMA_RESPONSE | INTERFERON_GAMMA_RESPONSE | 75 | 0.6768537 | 2.074181 | 2.107233e-05 | 0.0007492422 | 0.0004907317 | 551 | tags=51%, list=18%, signal=43% | Ifitm3/Ifitm2/Irf7/Bst2/Isg15/Pim1/Zbp1/Oasl1/Xaf1/Rnf213/Isg20/Eif2ak2/Nod1/Nfkbia/Herc6/Tdrd7/Trafd1/Mvp/Eif4e3/Parp14/Stat1/Socs3/Icam1/Nlrc5/Irf9/Epsti1/Myd88/Tnfaip3/Ifi35/Il10ra/Ifnar2/Psme2/Slamf7/Casp4/Ogfr/Znfx1/Casp1/Irf1 |
| TNFA_SIGNALING_VIA_NFKB | TNFA_SIGNALING_VIA_NFKB | TNFA_SIGNALING_VIA_NFKB | 91 | 0.6440568 | 1.992252 | 3.329965e-05 | 0.0007492422 | 0.0004907317 | 746 | tags=66%, list=24%, signal=51% | Tnf/Ccl4/Ptpre/Cebpb/Atf3/Klf10/Cd44/Ldlr/Sgk1/Sat1/Trib1/Dusp1/Egr1/Marcks/Nfkbia/Nr4a3/Egr2/Ier2/Fos/Traf1/Zfp36/Socs3/Pdlim5/Icam1/Klf4/Fosb/Rnf19b/Btg2/Dusp5/Tnfaip3/Jun/Nr4a2/Nr4a1/Maff/Mxd1/Pfkfb3/Nfe2l2/Kdm6b/Map2k3/Gadd45b/Map3k8/Atp2b1/Junb/Hes1/Irf1/Plek/Ptger4/B4galt5/Ppp1r15a/Ninj1/Sqstm1/Il6st/Mcl1/Tsc22d1/Smad3/Rhob/Stat5a/Dennd5a/Myc/Tank |
| INTERFERON_ALPHA_RESPONSE | INTERFERON_ALPHA_RESPONSE | INTERFERON_ALPHA_RESPONSE | 40 | 0.7355514 | 2.108775 | 9.441205e-05 | 0.0014161807 | 0.0009275570 | 342 | tags=40%, list=11%, signal=36% | Ifitm3/Ifitm2/Irf7/Bst2/Isg15/Oasl1/Isg20/Eif2ak2/Tent5a/Herc6/Tdrd7/Trafd1/Parp14/Irf9/Parp9/Epsti1 |
| IL6_JAK_STAT3_SIGNALING | IL6_JAK_STAT3_SIGNALING | IL6_JAK_STAT3_SIGNALING | 36 | 0.7292769 | 2.068579 | 1.512185e-04 | 0.0014367635 | 0.0009410381 | 375 | tags=47%, list=12%, signal=42% | Tnf/Il1r2/Csf2ra/Cd9/Pim1/Cd44/Cd36/Tnfrsf1b/Ifngr1/Ebi3/Il17ra/Pik3r5/Stat1/Socs3/Irf9/Myd88/Jun |
| APOPTOSIS | APOPTOSIS | APOPTOSIS | 50 | 0.6813605 | 2.005567 | 1.915685e-04 | 0.0014367635 | 0.0009410381 | 551 | tags=48%, list=18%, signal=40% | Ifitm3/Tnf/Lgals3/Atf3/Cd44/H1f0/Gpx1/Sat1/Isg20/Ifngr1/Pmaip1/Pea15a/Ank/Btg2/Jun/Dap/Sod1/Gadd45b/Tspo/Casp4/Ppt1/Rara/Casp1/Irf1 |
| COMPLEMENT | COMPLEMENT | COMPLEMENT | 58 | 0.6670681 | 1.983564 | 1.667369e-04 | 0.0014367635 | 0.0009410381 | 120 | tags=21%, list=4%, signal=20% | Fcer1g/Lgals3/Irf7/Ctsb/Msrb1/Cebpb/Anxa5/Pim1/Dgkh/Gngt2/Atox1/Cd36 |
In [67]:
path_out
'~/projects/2025/LeeNaeun/out/'
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [68]:
sessionInfo()
R version 4.4.2 (2024-10-31) Platform: x86_64-conda-linux-gnu Running under: Ubuntu 20.04.2 LTS Matrix products: default BLAS/LAPACK: /home/vincent/aprojects/Yeonjoo/anaconda3/envs/r_env_restored/lib/libopenblasp-r0.3.29.so; LAPACK version 3.12.0 locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=ko_KR.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=ko_KR.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=ko_KR.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=ko_KR.UTF-8 LC_IDENTIFICATION=C time zone: America/Los_Angeles tzcode source: system (glibc) attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] stringr_1.5.1 msigdbr_25.1.0 enrichplot_1.26.6 [4] clusterProfiler_4.14.6 scales_1.4.0 ggpointdensity_0.2.0 [7] future_1.58.0 harmony_1.2.3 Rcpp_1.1.0 [10] dplyr_1.1.4 ggplot2_3.5.2 Seurat_5.3.0 [13] SeuratObject_5.1.0 sp_2.2-0 loaded via a namespace (and not attached): [1] RcppAnnoy_0.0.22 splines_4.4.2 later_1.4.2 [4] pbdZMQ_0.3-14 ggplotify_0.1.2 tibble_3.3.0 [7] R.oo_1.27.1 polyclip_1.10-7 fastDummies_1.7.5 [10] lifecycle_1.0.4 globals_0.18.0 lattice_0.22-7 [13] MASS_7.3-65 magrittr_2.0.3 plotly_4.11.0 [16] ggtangle_0.0.7 httpuv_1.6.16 sctransform_0.4.2 [19] spam_2.11-1 spatstat.sparse_3.1-0 reticulate_1.42.0 [22] cowplot_1.1.3 pbapply_1.7-2 DBI_1.2.3 [25] RColorBrewer_1.1-3 abind_1.4-8 zlibbioc_1.52.0 [28] Rtsne_0.17 purrr_1.0.4 R.utils_2.13.0 [31] BiocGenerics_0.52.0 yulab.utils_0.2.0 GenomeInfoDbData_1.2.13 [34] IRanges_2.40.1 S4Vectors_0.44.0 ggrepel_0.9.6 [37] irlba_2.3.5.1 listenv_0.9.1 spatstat.utils_3.1-4 [40] tidytree_0.4.6 goftest_1.2-3 RSpectra_0.16-2 [43] spatstat.random_3.4-1 fitdistrplus_1.2-4 parallelly_1.45.0 [46] codetools_0.2-20 DOSE_4.0.1 tidyselect_1.2.1 [49] aplot_0.2.8 UCSC.utils_1.2.0 farver_2.1.2 [52] matrixStats_1.5.0 stats4_4.4.2 base64enc_0.1-3 [55] spatstat.explore_3.4-3 jsonlite_2.0.0 progressr_0.15.1 [58] ggridges_0.5.6 survival_3.8-3 systemfonts_1.2.3 [61] tools_4.4.2 treeio_1.30.0 ragg_1.3.3 [64] ica_1.0-3 glue_1.8.0 gridExtra_2.3 [67] qvalue_2.38.0 GenomeInfoDb_1.42.3 IRdisplay_1.1 [70] withr_3.0.2 fastmap_1.2.0 digest_0.6.37 [73] gridGraphics_0.5-1 R6_2.6.1 mime_0.13 [76] textshaping_1.0.1 scattermore_1.2 GO.db_3.20.0 [79] Cairo_1.6-2 tensor_1.5.1 spatstat.data_3.1-6 [82] RSQLite_2.4.1 R.methodsS3_1.8.2 RhpcBLASctl_0.23-42 [85] tidyr_1.3.1 generics_0.1.4 data.table_1.17.6 [88] httr_1.4.7 htmlwidgets_1.6.4 uwot_0.2.3 [91] pkgconfig_2.0.3 gtable_0.3.6 blob_1.2.4 [94] lmtest_0.9-40 XVector_0.46.0 htmltools_0.5.8.1 [97] fgsea_1.32.4 dotCall64_1.2 Biobase_2.66.0 [100] png_0.1-8 spatstat.univar_3.1-3 ggfun_0.1.9 [103] reshape2_1.4.4 uuid_1.2-1 curl_6.0.1 [106] nlme_3.1-168 repr_1.1.7 zoo_1.8-14 [109] cachem_1.1.0 KernSmooth_2.23-26 parallel_4.4.2 [112] miniUI_0.1.2 vipor_0.4.7 AnnotationDbi_1.68.0 [115] ggrastr_1.0.2 pillar_1.11.0 grid_4.4.2 [118] vctrs_0.6.5 RANN_2.6.2 promises_1.3.3 [121] xtable_1.8-4 cluster_2.1.8.1 beeswarm_0.4.0 [124] evaluate_1.0.4 cli_3.6.5 compiler_4.4.2 [127] rlang_1.1.6 crayon_1.5.3 future.apply_1.20.0 [130] labeling_0.4.3 plyr_1.8.9 fs_1.6.6 [133] ggbeeswarm_0.7.2 stringi_1.8.7 viridisLite_0.4.2 [136] deldir_2.0-4 BiocParallel_1.40.2 babelgene_22.9 [139] assertthat_0.2.1 Biostrings_2.74.1 lazyeval_0.2.2 [142] spatstat.geom_3.4-1 GOSemSim_2.32.0 Matrix_1.7-3 [145] IRkernel_1.3.2 RcppHNSW_0.6.0 patchwork_1.3.1 [148] bit64_4.6.0-1 KEGGREST_1.46.0 shiny_1.11.1 [151] ROCR_1.0-11 igraph_2.0.3 memoise_2.0.1 [154] ggtree_3.14.0 fastmatch_1.1-6 bit_4.6.0 [157] gson_0.1.0 ape_5.8-1
In [ ]: